LPTENS 03/30 



Field-theoretic approach to metastability in the contact process 

Christophe Deroulers^'0 and Remi Monasson^''^'0 

' Laboratoire de Physique Theorique de I'ENS, 24 rue Lhomond, 75231 Paris CEDEX 05, France 
^ CNRS-Laboratoire de Physique Theorique de I'ENS, 
24 rue Lhomond, 75231 Pans CEDEX 05, France 
CNRS-Laboratoire de Physique Theorique, 3 rue de I'Universite, 67000 Strasbourg, France 

(Dated: February 2, 2008) 

A 'quantum' field-theoretic formulation of the dynamics of the Contact Process on a regular graph 
of degree z is introduced. A perturbative calculation in powers of l/z of the effective potential for 
the density of particles (j){t) and an instantonic field ijiit) emerging from the formalism is performed. 
Corrections to the mean-field distribution of densities of particles in the out-of-equilibrium stationary 
state are derived in powers oil/ z. Results for typical (e.g. average density) and rare fluctuation (e.g. 
lifetime of the metastable state) properties are in very good agreement with numerical simulations 
carried out on D-dimensional hypercubic {z = 2D) and Cayley lattices. 

PACS numbers: 05.70.Ln, 02.50.-r, 03.50.Kk, 05.10.-a, 05.40.-a, 05.50.-|-q, 64.60.Cn, 64.60.My, 75.10.Hk, 
82.20.-W 



I. INTRODUCTION 
A. Motivations 

Recentyears have seen an upsurge of interest for the dynamical properties of out-of-equihbrium systems in statistical 
physics Systems of interacting elements are ubiquitous in physics and other fields, e.g. biology, computer science, 
economy... Most of the time the dynamical rules do not obey detailed balance or similar criteria which would ensure 
the existence of a well defined stationary distribution at large times. In other cases, a Gibbs measure does exist but 
is out-of-reach on experimental time scales, and all phenomena of interest e.g. the occurrence of dynamical phase 
transitions take place when the system is truly out-of-equilibrium. An example of such out-of-equilibrium phenomena, 
frequently encountered in condensed matter, in cellular automata or even in computationally-motivated problems is 
the occurrence of metastable states, or regions in phase space in which trapping may take place for a very long time 
before further evolution becomes possible. 

The calculation of the temporal properties of these systems often turns out to be very hard, even when dynamical 
rules look like innocuously simple. Over the past decade, however, various models and problems have been successfully 
investigated e.g. 0, 0, 0, IE 13 ■ Among the analytical methods used to tackle these systems, some rely on the 
observation that the master equation for a system of classical degrees of freedom may be seen as a Schrodinger 
equation (in imaginary time) where the quantum Hamiltonian encodes the evolution operator. Exact e.g. Bethe 
Ansatz Q or approximate e.g. variational or semi-classical techniques developed in the context of quantum field 
theory may be used to understand the dynamical properties of the original system 0- One important achievement 
made possible by this approach once combined with renormalization group techniques has been the calculation of 
decay exponents and the identification of universality classes in reaction-diffusion models [E |^ . 

The range of this "quantum" procedure is however not limited to the calculation of universal quantities. In this 
work, we show how it can be combined to diagrammatic techniques developed in the contexts of field theory and 
the statistical physics of disordered systems to quantitatively characterize the metastable properties of a well-known 
. !^ [ example of out-of-equilibrium system, the so-called contact process (CP) ^3 • In spite of its technicalities, this approach 
allows us to make predictions that can be successfully compared to numerical simulations. It is expected that it will 
$H ' permit to investigate metastability |llLll2j| . or other properties of various dynamical models. 
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FIG. 1: Profiles of tfie density p of particles versus time t for the Contact Process over a complete graph of A'^ vertices, initially 
filled with particles (p(0) = 1). A. Thermodynamic limit, — » cxd. From bottom to top: subcritical (A < Ac = 1, the density 
exponentially relaxes toward zero), critical (A — Ac, the density algebraically decays to zero as p{t) ~ t~^), and supercritical 
(A > Ac, the density exponentially relaxes to a finite value, p* = 1 — 1/A) cases. The density obeys a deterministic evolution 
equation I36II . and no fluctuation is present. B. Finite size lattice, with A'' = 100 sites. CP is a stochastic process, and the 
density profiles vary from run to run (we show here one run for each value of A). The density quickly relaxes to zero (subcritical 
regime) or a finite value (supercritical regime). In the latter case, the system is trapped in a metastable regime where the 
density fluctuates for a very long time around its plateau value (Inset, notice the difference of time scale), till a large fluctuation 
drives the density to the zero value. 



B. The Contact Process: Definition and Phenomenology 

We consider a regular graph G with vertex degree z and size N (number of vertices). Each vertex (or node, or site) 
may be empty or occupied by one particle. Hereafter, we focus on the continuous time version of CP where a particle 
is spontaneously annihilated with rate 1 independently from other sites, and an empty site becomes occupied with 
rate Xriocc/ z where riocc is the number of its occupied nearest neighbours |lCll |. 

The value of the parameter A strongly affects the behaviour of CP. For infinite size graphs e.g. infinite hypercubic 
lattices, there exists a critical value Ac of A such that 0,0], 

• if A < Ac, the number of particles (occupied sites) quickly decreases towards zero. Later, the system remains 
trapped in this empty configuration. 

• if A > Ac, the density p of particles reaches a plateau value p*(A) independently on the initial density |52l |. 



at criticality, that is, when A = Ac, the density eventually relaxes to zero with a slow algebraic decay p{t) ^ V 
This critical behaviour falls into the directed percolation universality class 0, ITgL ITtI I . Exponent a is equal to 
1 in dimensions larger than D,. = 4, and approximate expressions in powers of Dr — D in lower dimensions have 
been obtained through the use of renormalization group techniques 0, H, S Ull • 

These behaviors are displayed in Fig. The exact value of the critical parameter Ac is unknown in any dimension 
_D, but rigorous bounds and estimates have been derived p^ . 

For finite size graphs, the empty configuration, referred to as vacuum in the following, is an absorbing state for the 
dynamics. Starting from any initial configuration e.g. fully occupied state, CP will eventually end up in the vacuum 
configuration after a finite time tyac- This forces the above infinite size picture to be smeared out by fluctuations in 
the case of large but finite lattices mill El- Ac locates a cross-over between fast {tyac{N, X < Ac) — 0{log N)) 
and very slow {tyac{N^ A > Ac) expO(7V)) relaxations towards the vacuum configuration. In the latter regime, the 
plateau height, p*, merely defines an average value around which the density exhibits fluctuations until the system is 
driven to the vacuum through a very large fluctuation (Fig. ^i). On time scales 1 <C i <C tyadN, A > Ac), the system 
is trapped into a metastable state [20j. A (pseudo-)equilibrium probability measure for the density can be defined, 

P{p,N)^e^p{N7r*{p) + o{N)) . (1) 

Function tt*, which depends on A and other parameters e.g. the dimension D for hypercubic lattices, describes rare 
fluctuations from the average density p*. Its maximal value is zero for p = p* . Densities p distinct from the average 
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FIG. 2: The large deviation function n* for the density of particles in CP with parameter A = 2 for the D- dimensional 
hypercubic lattice. The D oo curve corresponds to the mean-field limit, and coincides with the case of the complete graph 
over N oo sites. Plot of the predicted value of vr*(p) at the order 1/D for D = 6, 3 and 2 are obtained with the expansion of 
Section IV. The nonconvexity of the curve for D = 2 shows the inaccuracy of the truncation to first order of our 1/D expansion 
in this range of densities. 

one are exponentially (in N) unlikely to be reached, and n* [p) < 0, see Fig. |5| In particular, the probability of a very 
large fluctuation annihilating all particles scales as exp(iV7r*(0)), and thus, 

W(iV) = exp(-Ar^*(0) + o(Ar)) . (2) 

The calculation of the large deviation function n* (p) is the main scope of this paper. For this purpose we use a path 
integral representation of tt* where particles are encoded into quantum hard core bosons, or 1/2-spins, and develop a 
diagrammatic self-consistent evaluation of the path integral which allows us to write a systematic expansion of tt* in 
powers of 1/z (Section II). In the infinite connectivity limit (z — > oo), this formalism reduces to the mean-field theory 
of CP, analyzed in Section III. Finite connectivity corrections to tt* are calculated in Section IV. As a by-product, 
we obtain an expansion for the critical parameter Xc{z) in powers of The validity of our calculation, effectively 
carried out up to order (and 1/z for tt*), is confirmed by numerical simulations performed on _D-dimensional 
hypercubic {z — 2D) and Cayley lattices. 



II. FIELD-THEORETIC FRAMEWORK 

A. Path integral formulation of the evolution operator 

We start by writing the master equation of CP using a quantum formalism, according to the familiar procedure of 
Felderhof , Doi and successors pH . |22 | . For each site i of the graph we define a hard-core boson with associated state 
vectors |0}i (empty) and \l)i (occupied), and creation and annihilation operators af and that anticommute on a 
single site but commute on different sites: [ai,af]+ = I, [a^,a^] = [ai,a~^] = [ai,aj] = (alternately, we could use 
spins 1/2 with a mere rewriting of the equations). To each occupation number = 0, 1 of site i is associated a vector 
\si). Then, to each state s of the graph (set (si, S2, . . . , s^v) of occupation numbers of all the sites) corresponds |23 a 
basis vector of a 2^-dimensional vector space, \s) = \si) ^ \s2) <^ ■ ■ ■ ® \sn) , and, to the time-dependent probability 
distribution P{s,t), the state vector \P{t)) = ^gP{s,t)\s). The master equation for P{s,t) is now equivalent to the 
evolution equation of the state vector 



±\P{t))=W\Pit)) 



(3) 
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where the evohition operator W is the infinitesimal generator of the transitions. For CP, W = Wai 

W^nn = J2il-a+)a, 

i 

i jei 

where j d i means that site j is one of the z nearest neighbours of site i. 

To map the stochastic process onto a path integral or field-theoretic formulation |^ |^ |^ , we introduce [25l 
continuously parametrized states suitable for hard-core bosons 53] . On each site of the graph, the state bra and ket 
are, respectively, 

{(t>,9\ = (1-0)^(01+05 exp(-z^) (1| , 

10,0) - (1-0)3 |0)-f 0^ eMiO) |1) (5) 
where G [0, 1], 6* G [0, 27r] and — —1. These states satisfy the closure relation: 

- /'d0 r"d0|0,0)(0,e| =1 . (6) 
Jo Jo 

To allow a simplification of the expressions in the translation table given below, we make use of 
-0 := — i ln[0/(l — (p)] + i 9 instead of 9, and introduce the following notations. 



(0,^/.| = (1-0)3 (^(O|+exp(-0.) (1| 

10,7/^) = (l-0)-3 (^(1-0) |O) + exp(^) |1) ) . (7) 



For the whole graph, a state is the tensor product of the states over all sites: |0,V') = '^iLi\4>ii'4'i)- Making use of 
Trotter formula |2q and of the closure identity ©, we obtain a path- integral expression for the matrix elements of 
the evolution operator exp(TM^) between times and T @,I3> 



4>(T)=4,tMT)=^t . . 

(0T,VT|exp(rW^)|0o,7Ao) = / P0(t)P#) exp -5[{0,7M] , (8) 

J ${o)=$oM^)=i'o ^ ' 

where the action reads 



C At If; 0.(0- 

•^0 U=l 



-5[{0,V}]^-/ At {Y.Ui)'^^'Wm)M))\ , (9) 



and the integral runs over all field configurations (f>{t),'ip(t) over the time interval t e [0,T] matching the required 
boundary conditions at initial and final times. 

Function W encodes the action of the evolution operator W on the states. Its expression is obtained by first writing 
W in normal order form thanks to the (anti)commutation relations, then using the translation Table U see 
For CP, we obtain W — Warm + AWcre with 

N 



Wan„(0,??) = ^0.(exp(V;,O - 1) , 

1=1 

1 ^ 

Wcrc(0,V') = -^0.^(1 -0,)(exp(-V;j)-l) . (10) 



z ■ 

1=1 jei 



The previous quantum formalism allows us to express the expectation value of any observable of interest. For 
instance, we may start at time t = from a random state \pq) with exactly A'o = po N occupied sites and project at 
the end on a state {pt\ with exactly Nt — pr N occupied sites: 

U):=T^E®-i((l-^»)|0) + -^»|l)) ' (PT|:=E®-if(l-^'')(0|+^''(l|) (11) 

\No) s ^ ^ s ^ ^ 
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operator in W 


expression in 


W 


1 


1 




a 


<jf> exp(?/;) 




a+ 


(1 - 0) exp(- 




a"*" a 








1-0 





TABLE L Translation table from operators in W into Lagrangian contribution to W. We have added a a"*" in the left column 
though this operator is not in normal order to show the consistency of the translation rules with the anticommutation relation. 



where the sums run over aU states s with Nq (resp. iV^) occupation numbers Si equal to 1, and the remaining N — Nq 
(resp. N — Nt) ones equal to 0. The probability that the density of particles equals px at time t — T given that is was 
equal to po at time t = is thus P{pT, T\po, 0) — {pt\ exp(T Vl^)|po)- Using the path integral formalism developed in 
this section, the logarithm 11 of this probability reads 



n(pT,r|po,o) = hi 



vmT^m exp (ptIWXT)) (0(O),V(O)|po) 



(12) 



where the boundary conditions for the fields at initial and final times are now free. Knowledge of the above function 
gives access to the large deviations function defined in Eq. ^ through 

7r*(p)= hm lim ln(p,T|p*,0) . (13) 

T^co N^Qo iV 

Notice that, although there is no notion of energy nor Hamiltonian in CP, the form of its path-integral formulation 
closely looks like a classical mechanics Lagrangian, with a kinetic energy term and an effective potential energy term. 

Calculation of the path integral on the r.h.s. of Eq. (|12p will be done through a perturbative expansion (in A) of the 
effective potential for the average values of the fields </>(<), ipit), following an approach used in the context of classical 
statistical mechanics ll^,!!^- This expansion allows us to calculate quantities of interest e.g. 0*(p), Xc, p*, ■■■ in 
powers of 1/D. In the following, we will closely follow the technique and notations of Ref. [33 which makes use of 
this perturbation expansion scheme to calculate equihbrium properties of the Ising model in large dimensions. The 
main difference (and complication) is that, here, fields depend on time. An application of this approach to the study 
of the dynamics of continuous spins models can be found in |3l{ . 



B. Diagrammatic expansion of the effective potential 

1. Constrained fields and conjugated sources. 

Let (t>i{t) and '0j(t) be two arbitrary functions depending on time t and site i, with 0j £ [0, 1], from which we define 
Xi{t) ■— (1 ~ (^„-(t))(exp(— ^,(t)) — 1). We choose as elementary (site-attached) operators 1 — (pi := I — afai and 
Xi := a^(I -f ai) — n |54|. and impose the constraints 



(I-0,(<)) = 1-0,(0 



(14) 



where (A) = {pT\exp{J^ dt'W) A exp{Jl^ dt' W)\po)/{pT\e^pijQ dt'W)\po) is the average value of operator A at 
time t. This can be done through the introduction of Lagrange multipliers (sources) in the evolution operator: W is 
changed to W'{t) + W"{t) 1 in the definition of {A} where 



E 



W"{t) := J2 



K{t) {l-c|>m-9^(t)X^{t) 



(15) 



Fields hi{t) and gi{t) are expected to be as regular as the imposed order parameters (t>i{t) and Xi(Oj ^ind are assumed 
to be (at least) once differentiable with continuous derivatives over the time interval t e]0; T[. However, to match with 
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the components of the final bra and initial ket, Dirac's ^-singularities may be present at t = and t = T. Note the 
introduction of a new parameter, /i, in front of the annihilation operator in the expression of W' . This parameter will 
result convenient for technical reasons only, and we will ultimately be interested in calculating quantities for /i = 1. 
This biased evolution operator allows us to express the logarithm of the probability that the final density equals px 
for a fixed set of order parameters, 

U[pT,T;{$,x}\po,0]^ln{pT\W{T,0)\po)+ [ dtW"{t) , (16) 

JQ 

and our task will be to compute VV(T, 0) := ex.p{J^ dtW'{t)). Requiring that be extremal with respect to (j>i(t) 
and (t) in addition to the constraints above ensures that hi (t) = gi (t) = at the extremum of 11. Therefore, at the 
saddle point, 11 in Eq. H16|l will coincide with 11 defined in Eq. I|12|) . 

The effective potential 11 can be expanded in a double power series in A, p, 

U= y2 p'Ua.b with Ua.b := ^{dxrid^)' n\x=^=o ■ (17) 

a,b>Q 

We calculate below IIo^o, that is, the effective potential in the absence of any evolution process albeit the one resulting 
from the kinetic constraint over the order parameters, and then expose how to obtain higher orders in a, b through 
a systematic diagrammatic expansion. A nice feature of this expansion scheme is that, at any given order a in A, we 
are able to resum the whole series in powers of p and, thus, to express our result as a unique power series, 

n = ^ A'' n,(/i) with n, := ^ / 0^,6 , (18) 

o>0 h>0 

and set p = 1 in the above expression. 



2. Calculation o/no,o. 

We set A = 0. W' decouples into a tensor product over the sites. The latter remain however coupled by the 
constraints that the bra {pt \ and the ket \po) correspond to configurations including exactly Nt — prN and TVq — poN 
particles respectively. We thus introduce two further Lagrange multipliers, vt and vqi to select the initial and final 
densities of particles. We replace {pt\ and \po) in Eq. (fT(^ with, respectively, {vt,Pt\ ■= (0|exp [i^Tj2iiPT — O't^i)] 
and ji/Q, Po) ■— {n) ^ sxp [j/q '^iio-t'^i ~ Po)] |0) where \0) — (|0) + |1))®^ is the sum of all possible configurations. 
Note that \i>q,pq) is normalized so as to represent a probability distribution. Once sites are decoupled, 11 may be 
expressed as a sum of site-dependent effective potentials, each depending upon 4>i{t)^ Xi(i), hi(t), gi{t), vt and vq. We 
will eventually optimize the resulting 11 over vt and to ensure that the final and initial densities are the requested 
ones. 

We send p to zero to make W' diagonal in the basis (|0),|1)). This allows us to compute exactly the evolu- 
tion operator VV(t2,ti) '■— exp(/j*^^ W'{t)dt)^ and then any average of operators or correlation function (CF) e.g. 

{a,{t2)al {ti)) {vt,Pt\W{TM) VV(^2,ii) a+ yV(ii, 0)|i/o, Po)- Evaluating ((I - ci)i){t)) and (x^), and imposing 
constraints (|14|l . we find back the rules listed in Table HI with over-barred fields i55||. The expressions for the sources 
hi{t) and gi{t) are then 56], for times < i < T, 

K{t) = ^^At) + (e^*'*) - 1) ^ln(l - 0,(t)) , 

g,{t) = Ain(l_0^(i)) , (19) 

from which we obtain hi{t){\ — 4>i{t)) — gi{t)Xi{t) = (1 ^ 0i(O)^V'i(i)- other words, the term in W" of Eq. l(TB|) 
gives back, aside boundary terms involving initial and final values, the 'kinetic' term of the action S (pj)l. 

It appears that the constraints (|14|l at times t — Q and t = T can be imposed by nonsingular sources hi{t) and 
gi{t) only if the required state vectors {4>i{T),tlji{T)\ and |0j(O), ^j(O)} are parallel to the initial bra {i/t,Pt\ and ket 
\votPo) respectively. To bypass this constraint, we introduce a singular term hifi5{t — 0) -|- hi^T5{t — T) in the source 
hi{t) — it is sufficient to modify hi(t) only and let gi{t) be regular, and we have also verified that singularities in hi(t) 
and giit) at times t g]0, r[ are absent unless (/'j(i) or V'i(*) are discontinuous, and that a discontinuity of the order 
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parameters is not favorable in terms of action and can be discarded. Optimization of IIo.o with respect to /i^.o and 
T yields 



/i,,o = V'»(0) + ln 



.(0) 



.(0) 



1^0 



(20) 



and allows to fulfill the constraints H14|) at initial and final times. 

Gathering all contributions to XIq^o and using Stirling's formula, we find after some algebra 



0,0 



E 



"^d^ ^^{t)^^ + Mpt - MT)) - ;^o(po - cj^M) - MO) In 0,(0) - (f - 0,(0)) ln(l - 0,(0)) 



+ N[po Inpo + (1 - Po) ln(l - Po)] 



(21) 



Notice that the sum of the last two terms in 0j (0) in Eq. H21|l is equal to the entropy of N noninteracting particles at 
density 0^(0). 



3. Perturbative expansions in powers of X and p. 

For the rest of this section, we call average of an operator A, and denote by (A), the ratio (vt, Pt\A W(r, 0)|i'o, Po) 
over (i^T, PT|W(r, 0)|fo, Po)- Let us introduce the operators 5*1 = — dt W^ann and 52 = — /g dt Were- These 
are directly related to 11 through the relations 9^11 = (— ^i) and d\Tl = {—S2), vahd for any A and p,. The 
average values of the operators and ^2 are {Si) = — /q dt X), 0i(*)(6xp('0i(i)) ~ 1) (for all A, 11) and {S2)\=o = 
— dt J2i 4>i (t) J2jiEi Xj (t) I (only for A = 0) respectively, and give the beginning of the expansion of 11. The sequel 
is obtained by iterative application of the following identities true for any (differentiable) operator A^ 



= (a^A> + (AC/), 5A(i) = (9AA) + (if) (22) 



with 

U = -Si 



(5i) + Tdt^ k/l,(i) {$^~Mt))+^^.9^{t){x^-X^{t))] (23) 
io , L J 

V = -52 + (52)+ / dtJ2 [dxh,{t) {4>, -Mt)) + ^x9^{t) {x^-X^{t))] (24) 



Though sites are coupled by V through 52, successive applications of operators U, V permit to evaluate the (derivatives 
of) CFs required to compute the expansion in powers of fi and A at A = where these CPs are factorized over sites. 
Notice that U, V and their derivatives with respect to fi and A that will appear in higher orders of perturbation 
involve the derivatives of the fields h and g. These can be expressed in a convenient way from the identities hi{t) — 
-d^^^^^n, giit) = -9^.(4)11, hi^T = -c^^(j.)n and /i,,o = "^^(o)^- In particular, dxhi^r = 9^^(y)(52) = and 
d\hiQ = (^.(o){S2) — 0, allowing us to remove all terms involving (5's in hi in the expressions of U and V. Also, 
d\hi{t) = c^.(t){52) and d\gi{t) = c^.(j)(52), so that the operator V writes in a convenient form when A vanishes: 
V\=o = Jq dtJ2i (i>iit)X]{t)/z, where 0i(<) = 0i - 0j(t) I and Xi(t) ^ Xi - Xiit) I express the deviations of the 
elementary operators with respect to their average values: {4'i{t)) = {xi{t)) = for all A,^. Finally, let us give some 
useful properties of U and V: 

1. (U) — (V) = (and if we write these operators in a natural way as sums over the graph sites, each summand 
{Ui) and (Vi) also vanishes); 

2. {uMt)) = {VMt)) = {UMt)) = (VUt)) = 0; 

3. Any CF of the type {Ai A2 ■ ■ ■ Ak U) vanishes if the time attached to the U operator is smaller or larger than all 
other times attached to the Ag operators. Indeed, a direct evaluation of (vtiPtIU and U\iyoiPa) after explicitly 
expressing [/ in a similar way as we did for Va=o shows that both vectors vanish. 

We now apply the above scheme to the calculation of 11. 
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III. MEAN-FIELD THEORY (z ^ oo LIMIT) 



In this Section, we first expose the analysis of the Contact Process on a complete graph. We then explain how this 
mean- field theory can be found back through the formalism developed in Section ITTI 



A. A simple derivation of mean-Held theory: the Contact Process on the complete graph 

Consider CP on the complete graph with N sites. As any two sites are adjacent, an exact account of the dynamics 
can be obtained from tracking the probability pn (t) that n sites are occupied at time t [s^ l . The master equation for 
these probabilities reads, 

^(t) = (n + l)p„+i(t) + ^^(n-l)(iV-n + l)p„_i(t)-(^n+^^n(iV-n)^p„(t) , (25) 

with the conventions P-i{t) — PN+iit) = 0. In the large N limit, we expect from the considerations of Section ITBl 
the following scaling behaviour for the probabilities [s^ . 

Pr^{t)=exp{NTTMF{n/N,t)) . (26) 

Inserting this scaling Ansatz into the master equation (|25|1 yields the equation of motion for the large deviation 
function TTMpip, t) of the density p of particles, 

dtTTMF (p, t) = Wmf [p, dpT^MF {p, t)] (27) 

where 

Wmf[</',V'] exp(^)-l)+ A (^(1-0) (exp(-^)-l) . (28) 

The very strong analogy between the above definition and the expression for the matrix elements of the creation and 
annihilation operators in Eq. IjlOl) will be explained in next Section. 

After a transient depending on the initial condition e.g. all sites are initially occupied and all densities but p = 1 
have zero probability, the large deviation function ttmf relaxes to its stationary value, 

^mf(p) = -^ + (1-p) (l-lnA-ln(l-p)) . (29) 

It is maximal and equal to zero in p = p*. The value in p = gives immediate access to the average lifetime of the 
metastable state with density p* , that is, the time it takes to the system to reach the empty state, see Eq. 



tvaci\N) ~ exp 



iV ( i + lnA - 1 



(30) 



This value may be successfully compared to the results of Fig. 4 in Ref. [321 . A direct numerical simulation of the 
evolution of the Pn{t) with N up to w 100 has allowed us to check the validity of scaling hypothesis (|26|) and to obtain 
a perfect agreement of the experimental distribution with Eq. (|29|l . 



B. Mean-field theory from the "quantum" formalism 

The first term, IIo^o, in the expansion of the effective potential 11 in powers of A and p is given by Eq. (|21|l . Then 
IIi^o = ^(5'i)o,o and IIo,! — —{82)0.0, the expressions of which in terms of (j> and x are given above. These are 
sufficient to establish the expressions of no(/i) and ni(/i), due to the third property of Section II. B. 3. Indeed, 

• Ila^o = for all a > 2. Proof: 9^11 — ~{Si U) since 9^5*1 = 0. The expression of (SiU) involves a double sum 
over the sites, say, of terms of the form {Ai Uj). Any such CF vanishes for all /i as stated in Section II. B. 3. 

• Ila.i = for all a > 1. Proof: df^dxH = —{S2 U) which, for A = and all p, reduces to a triple sum over the 
sites, say, i, j ^ i and k, of CF of the type {Ai 13 j Uk) where the times attached to Ai and Bj coincide. Again, 
these CF vanish due to the presence of the U operator. 
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As a result, setting /i = 1, we obtain 



Il[pT,T;{^,x}\po,0] = J2 



vt{pt - (t>i{T)) - vq{pq - (/'i(O)) 



0,(0) ln(0,(O)) - (1 - c^M) ln(l - 0.(0)) + Po Inpo + (1 - Po) ln(l - po) 

T 

dt 



+ I di U,(i)^(t) + 0.(i)(cxp(V.,(t))-l) + ^0,(t)5]x,(i) 



+ 0(A2) (31) 



For A = resp. for /i = 0, we are able to compute the action (|31|l directly, without using such a perturbative expansion 
as (|17|l . since the matrix of VF' is upper resp. lower triangular in the basis (|0), |1)). However, when both A and p are 
nonzero, we do not know how to proceed without the perturbative expansion. 

In search for a translationally invariant i. e. site-independent evolution, wc choose all quantities to be site indepen- 
dent and remove site indices. We also drop the bars over the fields to simplify notations. Equation 1)3 l|l then becomes 
ttmf = Tl/N with 

^MF [PT, T; {0, x}\po, 0] = ln(0(O)) - (1 - 0(0)) ln(l - 0(0)) + po Inpo + (1 - Po) ln(l - Po) 

+ i^TiPT - <i>{T)) ^ mpo ~ m) + dt(^m^{t) + WMF[m,m]^ (32) 

where Wmf is defined in Eq. (|28|l . The above expression suffices to the study the mean- field (MF) case i.e. when the 
site connectivity z goes to infinity. This may happen for large complete graphs, where each site is connected to all 
other sites {z = TV— 1), or in the D ^ oo limit of a Z?-dimensional regular lattice (z = 2D). As shown in Section IV. A, 
higher order terms in the A expansion give 0{l/z) additive contributions to 11 (within a site independent Ansatz), 
and vanish in the z — > oo limit. Therefore, ttmf is the exact mean-field expression for the action. To obtain the 
(exponentially in N) dominant trajectory of the order parameters (j){t),ip(t), we first extremize H32|) with respect to 
vttVoj to get the expected relations (j){T) = pT and 0(0) — pq. Notice that, if the initial state were a superposition 
of configurations with various densities e.g. dpo exp(iVg(po))|po), Eq- (|32|l would include an additive contribution 
qipo)', the most probable initial density, po, would then be given by the solution of dpqlpa) = ipiO). The resulting 
expression, 

r dt[m^(.t) + WMF[m,m]} , (33) 



TTMF = 

has then to be functionally extremized over the fields 0(t) and ip{t). The equations of motion (EM) for the fields are 
the Hamilton- Jacobi equations associated to Lagrangian 1)3311 . 

along with the already established initial and final conditions: 0(r) — px and 0(0) — pq. 

The solution of the EM yields the logarithm Nttmf of the probability to go from a configuration with density 0(0) 
at time to another configuration with density 0(r) at time T. To solve these equations, we make use of the fact 
that Wmf is a conserved quantity from Eq. (|34|l . which we denote by E. Then, the density and the time t can be 
expressed as functions of y := exp{ip), 

t= r dy'[{y'-Xf{y'-lf + 4EXy'{y'-l)]-'^^ . (35) 

The action of the trajectory equals the large deviation function, TTMF (p, T) — /J'(o)'^ dy' \n{y')dy(j>{y') + T E. The shape 
of the solutions of H34|l depends on the sign of "0(0). As shown in Fig.O if we exclude solutions where 0(f) is always 
zero jsj, three cases have to be distinguished: 

• if "0(0) = 0, the solutions have an infinite lifetime. tp{t) remains zero at all times [s^, and the field 0(i) obeys 
the first-order ordinary differential equation 

^W = -A0(O Ut)-^ + \) ■ (36) 
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This equation coincides with the mean- field equation for the density p{t), straightforwardly obtained when 
neglecting correlations between occupation numbers of neighboring sites. At large times, (j){t) tends to or 
p* = 1 — 1/A depending on whether A is smaller or larger than the critical value Ac = 1 as shown in Fig.Q'^- 
Notice that the action ttmf vanishes a.s ip — 0. 

• if "0(0) < 0, the solutions have a finite lifetime and 4>{t) 1 while ■^'{t) — oo (see Fig. |3J\&B). 

• if tp{0) > 0, the sign of the conserved quantity E matters. If > (which is necessary if A < 1), the lifetime 
is finite (see Fig. If £' < 0, the solutions are periodic (see Fig. |2f>). li E = 0, the situation is the natural 
intermediate between the two cases. In all cases, the lifetime (or the period) T{E) increases as \E\ diminishes. 



Nonzero fields ip select instantonic solutions allowing the system to escape the typical behaviour described by the 
— solution. These solutions have an extensive action, and are therefore exponentially suppressed when the 
system size increases. However, they are the solutions giving rise to large deviations of the density for a finite volume 
(Fig. n)3). The most probable fiuctuations correspond to long times solutions i.e. E ^ 0. In this case, T{E) diverges 
like logE', and the T E term vanishes. When A > 1, we have the simple identity (f){y) = 1 — y/A if < y < A, if 
y > X. In particular, the weight of the instantonic solution that goes from 0(0) = p* to (j){+oo) = p coincides with 
the stationary large deviation function for the density, 7r£{p(p), defined in H29|l . 

The large deviation function can be obtained at any finite time t too. To do so, we use expression (I32|l to express 
the probability of going from the state with distribution ^{p, t) at time t to the one the state with Tr{p, t + dt) at time 
t + dt. Special care must be paid to discretizing the term involving time derivations in a symmetric way between t and 
t -\- dt sls requested for path integrals '2^ . The resulting evolution equation for the large deviation function coincides 
with <(77|l as expected. 



C. Relationship with Martin-Siggia-Rose and Janssen-de Dominicis formahsms 

The above formalism is, to some extent, related to the treatment of Langevin equations by Martin, Siggia and Rose 
(MSR) and Janssen and de Dominicis jSSi^^ls^. If x{t) is a classical (scalar or vector) field, the Langevin 
equation 

^{t)^-V'{xit))+rjit) (37) 

describes its evolution in the potential V {V denotes dxV) under the random Gaussian force r]. rj is specified by 
its first two moments: ri{t) — 0, rj{t)ri{t') = 2TS{t — t') where T is the temperature and the overlines denote the 
noise average. Let P{xt,T\xo,0) be the probability that x equals xt at time T conditioned to its initial value. 
This probability can be expressed as a path integral through the introduction of a response variable x conjugated to 
X ,34. .35. .3&. ■37j ■ 



l.x(T)=XT r rt 

P{xT,T\xo,0) = Vx{t) Vx{t) exp / 

Jx(0)=Xn Jo 



fx(T)=x 
'x{0)=xo 

Functional optimization of H38(l with respect to x,x yields the classical equations of motion 



dT { ix{T) (^{t) + V'{x{t)) ] - Tx\t) 



(38) 



^^^V'{x{t))-^2Tx{t), ^ = mV'\x{t)) (39) 

from which we obtain that ^{t) — V'{x{t)) is a constant of motion. 

These EM are formally identical to those derived for CP H34I) in Section III {x playing the role of (p and x that of 
—itp) when << 1, with the following choice of potential and temperature, 

y'(0) = A0(0-r), 2r(0) = 0(i + A(i-0)) (40) 

Therefore, in the weakly fluctuating regime where ip is small and (p close to its most probable value p*, the system 
evolves in an effective potential V with fluctuations that can be described by a Langevin equation at temperature 
T |59l |. This description breaks down as soon as we consider large deviations, or transitions from the metastable to 
the empty state, as can be seen from the numerical comparison of the true solutions of Eq. H34|l with the solutions of 
Eq. 
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FIG. 3; Solutions of the mean-field equations of motion 1)34^ for A = 2 (p* = 1/2). A & B. For negative t/'(0) (= — 10~^ here), 
the solutions have a finite lifetime. After quickly reaching the most probable value of the density of particles, p* , cj>{t) stays 
for a long time in the neighborhood of p* but finally reaches unity, while 1/1(4) — » —00. The action ttmf tends to a negative 
finite value; the main contribution to it comes from the final jump of from p* to 1. C. For positive i/'(0) (= 10~* here) and 
E > Q, the lifetime is finite. After quickly reaching p*, <j){t) stays for a long time in its neighborhood but finally transits to the 
neighborhood of where it stays again for a long time, while xp{t) transits to In A and stays close to it. Finally ip{t) — > +00 
while cj>{t) 0. TTMF tends to a negative finite value, the main contribution of which is accumulated during the transit. D. 
For positive ^(0) > (= 10~* here) and E < 0, the solution is periodic. cf){t) and ip{t) oscillate between two values inside 
the [0, p*] and [0, ln(A)] intervals, respectively, and the action diverges to —00 by hops. E. Phase portrait with all types of 
solutions (oriented according to the time evolution). The phase space divides into four regions delimited by the solutions for 
which E = (represented in bold lines); in each region the sign of E is indicated. The quasistationary state appears on this 
diagram as the crossing point (tf) — 0, tf) = p*) of two solutions with E = 0: it is stable along the <j!> direction but instable along 
the ip direction. 
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Extending the validity of the Langevin equation approach to the full domain of <j) and would require the use of 
a MSR-like formalism [3 HE HE US • Interestingly, MSR stated in their original work that the knowledge of the 
physical field (j) alone is not sufficient to compute all quantities of interest beyond the Gaussian approximation, and 
that the introduction of a second operator (corresponding to our "0, or to a and , whereas cj) corresponds to a^a) 
to express the response functions of the physical field was required. 



IV. FINITE-DIMENSION THEORY {l/z EXPANSION) 



A. Analytical calculation 



Calculation of the corrections to mean-field theory requires the knowledge of higher-order terms in the expansion of 
n in powers of A and ^. Strictly speaking, our expansion is, after we resum the /i-expansion, an expansion in powers 
of A, or, more precisely, of A/z. It naturally gives access to an expansion of 11 in powers of l/z as shown below. 



1. The diagrammatic expansion 



To show how this expansion works in practice, let us consider the first correction to Eq. (|32|l . 2112,0 — &xI\\ofl — 
~d\{S2)ofi = —{S2V)ofi = (V^^)o,o since 8x82 — and {S2V) = —{V"^). Vo,o involves a sum over couples of 
neighboring sites that can be written as a sum over oriented links. The term Xj will be represented by a link going 
from site i to site j. Therefore (F^)o,o writes as a sum over couples of oriented links. When A = 0, sites decouple and 
terms where a site carries a single operator cj) or x vanish (Section II. B. 3). We are thus left with terms where two 
parallel links loop over two neighboring sites i,j on the lattice, 

«<t>J= / dt [ dt' {Mt)xAt)Mt')xj{t'))o,o , (41) 



and terms with the same structure but links pointing in opposite directions. As A = 0, this four-point CF factorizes 
into a product of two two-point CFs: {(f>i{t) Xj{t) (j)i{t') Xj{t'))o.o — {<l>i{t) ^i(i'))o,o (Xj(^) Xj(^'))o,o- Each of these CFs 
can be straightforwardly computed e.g. (0i(t) 0i(i'))o,o = (1 — 0j(^2)) where ti = min(t, t') and t2 = max(<, t'). 

The final result reads 

^ <5>J = -2 / dt2 dt, 0,(ii) (1-^,(^2)) (1 + X,(ii)) X,it2) ■ (42) 
Jo Jo 

The antiparallel two-site loop diagram may be evaluated in the same way. Notice that the final outcome does not 
depend on the indices i,j of the sites, provided these are neighbors on the lattice. We will therefore drop site indices 
in the following. Then, we count the multiplicity of each diagram, equal to A'^ z for both parallel and antiparallel 
two-sites loops. Finally, each diagram gets a factor 1/z^ coming from the (A/z)^ factor in Were- As a result, the net 
contribution will be of the order of l/z, that is, l/I? on a /^-dimensional hypercubic lattice. 

We are now able to write the general expression of the l/z expansion in a graphical way. For instance, for a 
hypercubic lattice of dimension D (and site connectivity z = 2D), 

" = '"'+2!(2D) 



-D^+... (43) 



3!(2i:i)3 4!(2i:>)4 ' ' A\{2Dy 

A^ 

Each undirected diagram in the above expansion represents the sum of the correlation functions (like the one entering 
112.0 Sind evaluated in the previous paragraph) that share the same support on the graph but differ by the orientations 
of the links. There are 2^ link orientations for an undirected diagram with £ links. The coefficients in front of the 
diagrams take into account the power of A/z (equal to the number of links in the diagram), the inverse factorial from 
the Taylor formula, and a combinatorial factor. This combinatorial coefficient is the product of the multiplicity of the 
undirected diagram (number of times it can be drawn on the lattice, divided by the lattice size N) and of the number 
of ways to associate a time t to each link of the diagram when evaluating the CF |60| . As one may infer from the first 
diagrams, there is only a finite number of terms contributing to a given order in l/z. 
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2. Summation of the fi-expansion and memory kernels 



n being the logarithm of a generating function of the fields, all the diagrams entering expansion H43|l are connected. 
Moreover, for = all nonirreducible diagrams i.e. which may be cut into two, or more, pieces by removal of a vertex 
vanish, as in the virial expansion (3^ and possibly the Ising model |30j | . This statement does unfortunately not extend 
to nonzero /i. To get the /x-expansion we make repeated uses of the leftmost identity in Eq. (|22|l . As df^U — d^j^V = 0, 
this amounts to insert U operators (one for each power of /x) in the CFs, yielding contributions of the form 

Dn .^ I dti [ dt2... [ dtn+2{Mtl)Mt2)Mt3)---Mtn+l) Mtn+2)) ■ (44) 



"'0 



As stated in Section II. B. 3, such a term vanishes if one of the times t2, ts, . . .tn+i associated to the U operators is 
the minimum or maximum of all times. If not, the integrand is the product of factors 

- expV>(tfc) 

for each lJ{tk) and a factor depending on the other operators involved in the CF. If there are strictly less than four 
other operators in the CF, or no U operator between the second and the second last times, this latter factor is equal 
to the CF where all U have been removed; for instance, the term £)„ defined in Eq. (|44|l reads 

Dn^2 dtj^ dt' (Mt) Ut')) (^^^ dt" W)^ . (46) 

If the original CF involves more operators, or a 'badly' located [/, there appears a kind of disentanglement of 
the original operators. Consider for instance the four-field CF D' :— {(t)i{ti)(l>i{t2)<pi{t3)(j)i{t4)) with ti < <2 < 
^3 < t4. A direct evaluation shows that D' equals — 0j(ii)(— 1 + (/!'j(t4))(— (/)j(t3) + 1 + 3(j)i{t2)4>i{t3) — 2cf)^{t2)) and 
cannot be expressed as a product of four factors, each of which would be associated to a time tj. We shall say 
that D' is entangled, here due to the f2,i3 term. Insertion of at least one Ui between t2 and ta, as in D" := 

{Utl)Ht'l)Mt'2) ■ ■ ■ Mt'nJMh)Mtn, + l)Uiit'n,+2) ■ ■ ■ Mt'nJMt3)Mt'n, + l)Mt'n,+2) ■ ■ ■ MtnJMU)) with h < 

t[ < . . . < i'j^ < t2 < t'ni+1 < . ■ ■ < t'n2 < ts < t'n2+i < ■ ■ ■ < Kia < ^4 rcmovcs this entanglement. D" is the product 
of factors ^i(t'j.), one for each U{t'j.), and of the disentangled (factorized) expression — 0j(ti)(— 1 + (/)j(t4))(— 1 + 
2^i(t2))(-l + 2^,(^3)). We conclude that a CF of the type 

pT pT pT pT pT pT pT 

D',[:^ dU dU dt2 dti dt[ d^.../ d4(0,(ti)^,(t2)0,(i3)^.(t4);/»(<i);/»(i2)---?^.(4)) (47) 



3 f dt2 1 dti (^{Mti)Mt2)Mt3WM)) - {Mti)Mt2)Mt3)Mt4))dis) 



Jo Jo Jo Jo Jo Jo Jo 
can be expressed as 

D'; = 24 / dti f dts [ d<2 / dti{Mti)Mt2)Mh)MU))dis( [ dt'^,{t) 
Jo Jo Jo Jo 

pT pt4 pt^ pt2 

+ 24 / dti I dta / dt2 

Jo Jo Jo Jo 

rt2 _ Z"*" _ \ " 

dtUt)+ dtUt)) (48) 

tl Jt3 / 

where the 'dis' subscript means that the disentangled expression for the CF integrand (. . .) must be taken. 

It appears that, at any given order in A, the series expansion in powers of fi can be easily summed up, yielding 
memory kernels in the resulting liaifJ-) with a > 2. For instance, Jq dt Jq dt' {^i{t)4>i{t')) is replaced with 

J2f^D,,= I dt I dt' mmt')) exp ( / dt" Un 1 , (49) 

n>2 



Jo \ Juiin{tJ' 



the memory kernel being the exponential term, where fi is eventually set to one. Due to the presence of the kernel, 
the CF between two operators at times t and t' decreases with a time difference |t — i'| (remember ^ < 0). 
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This property extends to more than two operators. The contribution of a diagram with n Unks, thus of order A" in 
the A-expansion of H, may be written as the sum of a finite number of terms of the form 

F:^ ( dt„ /" " dtn-i ... [' dhP exp III f ' dt^it) + . . . + i„ /" " dtlit)) (50) 
^0 Jo Jo \ Jtl Jt„-i J 

where P is a polynomial in (/'(ti), exp(— . . . , and exp(— ■!/)(t„)) and ii, 12, . . . , in are positive integers. 

The presence of the memory kernels ensures that the integrand sharply decreases as <„ — ti increases. 

One may wonder why memory kernels appear in the expressions above, whereas CP is a Markovian stochastic 
process. This is the consequence of the projection of the complete distribution of states, \P{t)), that is, the knowledge 
of the probability of occupation or vacancy of all sites, onto a partial description where we keep track of the order 
parameters 4'{t)^ '0(t) only. Beyond mean field, the degrees of freedom which were discarded pop up as non-Markovian 
contributions to the dynamical evolution of the order parameters. This phenomenon is well known, and can be 
illustrated with an elementary example proposed in Appendix A. 



B. Results and comparison with simulations 

Following the above recipes, the effective potential tt may be expanded in powers of 1/z, see Eq. (|43|) . The explicit 
expression for the first correction to mean field, coming from the two-site loop diagram, reads 

TTi = £ dt^il - q^{t2)f (e-^(*^)-l)^*'dti</>(ti) ((l-<^(ti))e-'^(*^)+<^(ti)) exp (^2 dt"e(t")) (51) 

The 0(l/z^) corrections to tt required the calculation of the diagrams listed in Eq. (|43() . The resulting expression, 
TT2, for the D-dimensional hypercubic lattice (where z = 21?) is too lengthy to be given here, but can be obtained 
with the help of a computer algebra software. This gives, as a by-product, the expression 7r2 for the Cayley lattice 
(infinite graph without loops where all sites have exactly z neighbours) after removal of the square diagram luj.. 

Functional optimization of the resulting expression for tt with respect to (t){t) and ?/'(t) yields EM which include 
corrections of orders 1/z and 1/z^ to the Hamilton- Jacobi equations H34|) . These corrections involve terms with 
multiple integrals on the time. 



1. Corrections to the density p* of particles and the critical parameter Ac 

We first concentrate on the solution to the EM with vanishing tp(t) at all times t. The equation for (j){t) then 
simplifies; for instance, to the first order in 1/z, we obtain 



dt 



{t) = -XHt) (^0(t)-l + l^-^(l-0(i))^ ^*dt'0(O exp(^-2^*^^^) . (52) 



With the help of computer algebra, the EM for (p{t) can be written up to order 1/z^, and its asymptotic behaviour 
analyzed. We find that 0(t 00) equals or p* > depending on the value of the parameter A with respect to its 
critical value Ac. The asymptotic density reads 

* . 1 1 6A2-H11A-I-3 , 3, 

^^1-A-A^ 6A^ + ^(^/^) ^''^ 

on the D-dimensional hypercubic lattice (z = 2D), and 

* . 1 1 6A2 + llA-3 , 3, 

^^1-A-A^ 6A^ + ^(^/^) ^''^ 

on the Cayley tree where all sites have z neighbours. The critical value Ac can be easily obtained from the above 
equations as the value of the parameter A below which no state with a finite density of particles can survive. This 
is intended to be the lower critical value Aci„„„, on the Cayley tree 0|, where there exists an intermediate range 
[•^Glower ' '^Cuppor ] valucs of A where nonuniform metastable states can survive without invading the whole graph 
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FIG. 4: Numerical data for the second-order correction in Xjz to Ac (2) for the Cayley tree (left panel) and the hypercubic 
lattice (right panel). The critical value Ac is obtained from numerical experiments on the conserved contact process (GPP), 
see main text. Data are plotted as 7? (Ac(z) — 1 — 1/z) versus I/2. The intercept a with the vertical axis (z — > 00) compare 
very well with the theoretical predictions a = 4/3 — Eq. 15611 . left — and 7/3 — Eq. 15511 . right. Dashed lines are tentative 
linear fits (with fixed origin at 1/z = 0) of the data, whose slopes should be given by the 0(1/2^) expansion. The estimate of A 
for each z was obtained through an extrapolation of data for finite numbers of particles A'^ (up to 4000) to infinite A''. For each 
value of A'^, we ran 10000 simulations and estimated error bars from the statistical fluctuations. The upper and lower values of 
Ac(2; A'^ = 00) were then obtained through a linear fit of Ac(2; 1/A'), and plotted at the extremities of the error bars. 



while, above Acupp„,. , there is a single uniform metastable state — on hypercubic lattices, both thresholds coincide. 
Setting p* = 0, we obtain 

Ac - 1 + - + + 0(1/^') (55) 

on a the hypercubic lattice, and 

Ac = 1 + - + + 0(1/^^) (56) 

Z 6Z^ 

on the Cayley tree. These results are compatible with the rigorous bounds (1 — 1/{2D))~^ < Ac < 4 (hypercubic 
lattice) and 1 < Ac < (1 - 2/z)-^ (on the Cayley tree) [lll3|- Notice that the lower bound Ac > (1 - 1/{2D))-^ 
was obtained |T^ through a two-site calculation. When discarding all diagrams but the two-site loops, we find 
Ac = 1 + 1/{2D) + 1/{2DY + 0(1/(21))^), that is, the same bound. For small dimensions, our asymptotic results are 
of poor quality. We refer the interested reader to the various works that give more precise estimates for Ac e.g. for 
hypercubic lattices in 1 < D < 5 ref. [s^ , for D — 1 through series expansions ref. |40l |4]| . 

Our values of p* are in good agreement with numerical simulations carried out on hypercubic lattices in dimensions 
up to D = 10 for several values of A (= 1.5,2,3). Simulating large size lattices in higher dimensions would require 
prohibitive memory space. Instead, we have performed simulations of the conserved contact process (CCP) 42] . This 
is a canonical counterpart of CP where the number of occupied sites is kept constant, but A fluctuates. The stationary 
properties of both processes are, in the limit of infinite lattice size, equivalent |4j,|4j. In particular, simulating CCP 
with a (not too large) number N of particles on an almost infinite lattice (in our simulation, of size L — 2^^) amounts 
to simulate CP with a vanishing density. The average value of A in CCP then coincides with the critical Ac in CP. We 
have been able to simulate CCP on hypercubic lattices with up to = 4000 sites and in dimensions up to £> = 80, 
or up to z = 81 on the Cayley tree |61| . Results are displayed in Fig. ^ and are in very good agreement with our 
theoretical predictions for Ac (I55I56|I and previous simulations |39||. 



16 



2. Corrections to the distribution of particle densities 



The analytical resolution of the EM when ^{t) ^ is difRcult. We have restricted ourselves to a first order in 1/z 
expansion around the infinite lifetime solutions for which we have an analytical expression in the mean-field case. It 
is convenient to parametrize the real time t in terms of yo introduced in Eq. (|35|l : j/o = 6xp('(/'o(i)) where ipo{t) is the 
mean- field expression for ipit)- In thermodynamic limit N oo with E (infinite lifetime), the t —>■ +oo limit 
translates into the yo X limit. 

In this setting, we may write (j){yo) = 0o(yo)+0i(2/o)/^ + O(l/z^) and similarly ip{yo) = ■0o(2/o)+^i(yo)/-z+O(l/z^), 
where <j)o and tpo were computed previously for the mean-field case. From the EM, we derive two coupled first-order 
linear differential equations for 0i and ipi: 



(,„ - - JL U 1 -(» - - ; + 1 - »)(» + 1) *.(») I + (57) 

dyo\Myo) yo \ 2A(yo-l) (yo-l)^ + A-l \ Myo) 
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/ (60, 

where a is defined by A =: 1 + 2/q;. The solutions of these equations diverge in yo — I for all initial conditions on 
4>i,tpi in ?/o = 1 except for <j>i{l) — — 1/(2A^), "01(1) — which precisely amounts to setting ■0(2/0 = 1) = and 
4>{yo = 1) = p* up to order 1/z^. We choose these initial conditions in the following. 

The resolution of the equations for (pi and ipi can be done in part analytically and in part numerically. First, we 
treat the neighborhood of yo to characterize exactly the (nondivergent) singularity in this point. We have calculated 
the solutions up to the order {yo — l)^ In |j/o — 1| included, which is sufficient to obtain the values of (/)i(l -I- e), '01(1 -I- e) 
slightly off (below or above) the singularity (with e = ±0.001 to 0.003) with a good numerical accuracy. Such an 
expansion amounts to a short-time expansion if one starts ai t ~ from the typical state where the density of full 
sites is p* and y ~ 1. Then, we start to solve the differential equations from the value yo ^ 1 + e using a Runge-Kutta- 
Fehlberg procedure. Note that the coefficients of the EM for 0i and tpi are rather simple, but the nonvanishing second 
member involves, for generic values of A, hypergeometric functions. To simplify the numerical resolution and to make 
it more precise, we restricted ourselves to the case A = 1 + 2/a with a a positive integer, where these hypergeometric 
functions reduce to polynomials and logarithms. 

This pcrturbative resolution yields a parametric representation of the large deviation distribution, tt* (p) — 7rJ[p (p) + 
7tI{p)/z + 0{\/ z^) with parameter yo- More precisely, we have 



p(yo) (P{yo)^Myo) + ^^^^ + o{\/z^ 

niyo) = plnyo-Soiyo)- — Si{yo) + 0{l/z^) (61) 



z 



where Soiyo) = Int/o - (z/o - 1)/A and 

X^ Ji X-yi\yi-lJ J I y2 - 1 \X-y2J 



The resulting curves for 7T*{p) are presented in Fig. [3 Apart from TT*{p*) = given by 1)531 154|l and reached for 
2/0 = 1, another point of interest may be located analytically, namely 7r*(p — 0), reached for yo = X. This is related 
to the lifetime of the metastable state, t^ac ~ exp{— Ntt* {p — 0)). Some values are listed in Table Hll for integer a. 

To check the accuracy of our perturbative expansion for tt* , we have performed simulations of CP on 6-dimensional 
hypercubic lattices with periodic boundary conditions. Sizes N = with L ranging from 3 to 6 are large enough 
that the system gets trapped for a time greater than the simulation run into the metastable state. We simulated the 
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A 


7r(p = 0) 


Decimal approximation 


3 


- In 3 + 2/3 + (-53/18 + ti'^/S) /z 


-0.432 - 


h 0.345/2 


2 


- In 2 + 1/2 + (-107/36 + 7rV3)/2 


-0.193- 


h0.318/z 


5/3 


- ln(5/3) + 2/5 + (-5413/1800 + ■n'^/'i)/z 


-0.111 - 


h 0.283/^ 


3/2 


- ln(3/2) + 1/3 + (-1823/600 + ir'^/S) /z 


-0.072 - 


h 0.252/2 


7/5 


- ln(7/5) + 2/7 + (-270281/88200 + 7rV3)/z 


-0.051 - 


h 0.225/^ 



TABLE 11: Some analytical values for the 1/z expansion of 7r(p = 0), the logarithm of the probability of reaching a configuration 
with vanishing density in the metastable state of the CP on a regular graph of coordination number z. Analytical calculation 
can be performed for values of A = 1 + 2/a with a a positive integer only, but 7r(p = 0) can be numerically estimated to order 
1/z for other values of A. 




P 



FIG. 5: Comparison, for D = 6 and A = 2, between the 1/73 predictions for the large deviation function 7r*(/9) and numerical 
results for CP on hypercubic lattices with periodic boundary conditions of linear sizes L — 3, 4 and 6. All curves were 
horizontally and vertically translated so that their maxima have the same coordinates p*,n* = 0. Due to the relatively high 
dimensionality of these systems, the range of values of p explored during the simulations is very concentrated around p* unless 
L is small. Note that the curves for different lattice sizes, once translated, seem to depend only weakly on L. Inset: blow-up of 
the top region; data compare well with the 1/D theoretical result (continuous line), and are clearly distinct from the mean-field 
result (dashed-dotted line). 

system starting from p <C p* or p ^ p* until it reached equilibrium. The equilibration times (expressed as the number 
of elementary steps) was found to correspond to continuous times .451] t^q ranging from 25 to 60. Then we run again 
the simulation for times t — M t^q with very large values of Af going from AI — 2.5 x 10^ for L = 6 to M = 5 x 10^ 
for L — 3, and recorded the histogram of p over this time interval. This gives a very good approximation of the 
quasistationary distribution n* since the system was already equilibrated. Numerical results for 7r*(p) are presented 
in Fig. |S| The maximum of it* (p) vanishes in the thermodynamic limit only, and the value p* at which this maximum 
is reached is also subject to finite-size corrections. To make the comparison easier, we have vertically and horizontally 
shifted the experimental curves so that they all reach zero at the same value of p. Once this translation was done, 
the numerical and theoretical curves are in very good agreement, and are easily distinguishable from the mean- field 
curve (Fig.lSJ. This shows that the curvature (unaffected by the translations) of 7r*(p) and, hence, the amplitude of 
the fluctuations in the metastable state and its lifetime, is correctly predicted by our 1/z calculation. 
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V. CONCLUSION 



In this paper, we have calculated the out-of-equilibrium distribution of densities of particles in the Contact Process 
on large regular lattices with degree z using a quantum field-theoretic formulation for the evolution operator. The 
calculation is based on a perturbative calculation of the effective potential for the density (\)(t') and an instantonic 
field ■(/'(i) as function of time t. Interestingly, this instantonic field naturally emerges from the parametrization of the 
quantum hard bosonic (or spins-1/2) states needed to represent the sets of occupation numbers. 

The finite connectivity corrections to the mean-field case {z = oo) lead to the appearance of memory kernels 
compensating the loss of information due to the tracking with time t of global fields ^'(0 only. Though calculations 
may rapidly become involved, we have been able to show the very good agreement of the predictions they give with 
numerical experiments on the average density of particles and, more generally, the whole distribution of densities in 
the out-of-equilibrium metastable state of CP. 

While the use of a quantum formalism was known to be very efficient to access universal quantities e.g. the decay 
exponent of the density with time at criticality, the present study shows that nonuniversal quantities can be calculated 
too. We hope that this approach will turn out to be useful to the analysis of the numerous far-from-equilibrium 
systems encountered in physics or related fields. From this point of view, a remote but promising field of applications 
could be the analysis of algorithms in computer science where out-of-equilibrium dynamics over (discrete) variables 
abound, and metastability phenomena are present |46l |4^ . Hopefully our approach will also permit to complete 
the average-case analysis of backtracking algorithms initiated in through the systematic control of the non- 

Markovian effects ignored so far ^] . 
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APPENDIX A: EMERGENCE OF MEMORY KERNEL WITH HIDDEN DEGREES OF FREEDOM 

Consider two variables x(t), y{t) obeying the Markovian evolution equations 

J(<) ^ -x{t)+(iy{t) , (Al) 

^(t) ^ -ay{t)+x{t) , (A2) 

with initial conditions a;(0) = 1, ?;(0) ~ 0. This system can be easily solved to give x and y as functions of time t. 
Assume instead we want to write an evolution equation for x only. Solving Eq. (|A2|) . and plugging the resulting y{t) 
in Eq. I|A1|I . we obtain 

J(i) = -xit) + /3 ^* § e-"(*-*') x{t') . (A3) 

As a result of the existence of a hidden degree of freedom y, the effective equation on x is not Markovian when 0, 
and includes a memory kernel whose time constant is that of the y variable. 
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